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We study propagating mode-I fracture in two dimensional amorphous materials using atomistic 
simulations. We used the continuous random network (CRN) model of an amorphous material, 
creating samples using a two dimensional analogue of the WWW (Wooten, Winer & Weaire) 
Monte-Carlo algorithm. For modeling fracture, molecular-dynamics simulations were run on the 
resulting samples. The results of our simulations reproduce the main experimental features. In 
addition to achieving a steady-state crack under a constant driving displacement (which had not yet 
been achieved by other atomistic models for amorphous materials), the runs show micro-branching, 
which increases with driving, transitioning to macro-branching for the largest drivings. Beside the 
qualitative visual similarity of the simulated cracks to experiment, the simulation also succeeds in 
explaining the experimentally observed oscillations of the crack velocity. 

PACS numbers: 62.20.mm, 46.50.+a 



The problem of dynamic fracture has been of great interest over the last two decades. Experiments, most of them 
done in amorphous materials [l|-[sj, have shown the existence of steady-state propagating mode-I fracture, under an 
imposed constant driving strain. The most interesting results concern the critical velocity at which a steady-state 
crack stops propagating on the sample midline, and starts either to bifurcate or to create micro-branching. The 
onset of this behavior was found to occur at velocities dramatically smaller than those predicted by Linear Elasticity 
Fracture Mechanics 0,|8|- 

Extensive theoretical efforts have been devoted to a variety of models in order to understand the high velocity regime 
instability. In lattice materials, steady- state cracks and the onset of instability have been observed in simulations 
and the corresponding solution for the steady-state crack [ol, [lo[ has been obtained, When the steady-state solution 
becomes linearly unstable, the crack starts to oscillate or to bifurcate . However, the behavior past the instability 
point did not correspond closely to the experiments [HI, |l5| . This is presumably due to the mismatch between the 
anisotropic lattice employed in the models and the disordered isotropic (on long scales) nature of the experimental 
amorphous materials. There has been some success in exploring the instability using phenomenological mesoscale 



171 , Il8j ; these however are difficult to quantitatively 



approaches based on the phase-field [16| or on cohesive zones 
relate to an underlying microscopic picture. 

Here, we focus instead on isotropic atomistic models. There have been a few earlier attempts in this direction [lol - 
[22I . These have focused on one standard model for glass, the binary-alloy model] here, particles of two sizes with a 
2-body (often Lenard- Jones) potential are mixed so as to avoid crystallization. Unfortunately, these kind of models 
do not allow for freely propagating cracks at all, under a constant driving. Only by continuously increasing the 
driving does the crack propagate, in direct contrast to experiment. In retrospect, this is not a surprising result, as 
this class of models yield large plastic zones. Thus, these mixtures can be used for metallic glasses, but not for brittle 
materials. 

Instead, we turn here to a different kind of model for describing amorphous materials, the "continuous random 
network" (CRN). This model goes back to Zachariasen, [23| who proposed it as a model for glass. In the CRN, the 
amorphous material is described as a "corrupted" ordered material. Wooten, Winer and Weaire 23, [25| (WWW) 



presented an algorithm for generating high-quality CRN's using Monte-Carlo techniques. Much effort has been 
expended in generating large 3D CRN's with thermodynamic properties matching those of several specific amorphous 
materials [26]. We note that Lu et. al. [22] has studied fracture in a small scale (1854 atoms) CRN model. However, 
they generated their network from a Poisson distribution of mass points, resulting in large voids, and their model did 
not exhibit steady-state fracture, but rather void coalescence. 

We generated two-dimensional CRN's by a 2D-analogue of the WWW algorithm (with 162 x 272 ^ 44, 000 atoms 
in the largest sample). The potential used in this constuction included both a 2-body central force and a 3-body 
bond-bending force: |27j : 
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FIG. 1: (a) An image of a steady- state crack in the CRN model. We see (in the circle) the origin of micro-branching. The 
crack reaches a "tough zone" and breaks parallel bonds, but it does not have enough energy to create a visible micro-branch. 
?7 = 2, ^/ag = 2.7 (b) At higher driving (77 = 1.5, ^/ag = 3.1), when the crack reaches the same tough zone, (displaced in this 
figure to enable showing its downstream development) the crack now has enough energy to create two branches. Eventually 
one branch wins and continues as the main crack. 



with Oc = 120° and a = 4. In our simulations, every atom has exactly 3 neighbors. 

Once we have constructed an appropriate amorphous material, we use an Euler scheme molecular-dynamics simu- 
lation for the fracture study. We use a piecewise-linear force law, in which the force decreases to zero immediately 
when the radial distance between two atoms exceeds a certain threshold (which we set at £ = 1.25a). The variation of 
bond length in the material is sufficiently small that no bonds are broken in the uniformly stressed state, before the 
initial crack is seeded. In addition, we used a Kelvin-type viscosity (with a viscose parameter rj) [13, 28, 29], and the 
3-body force law - the gradient of Eq. (TJ that is removed when any of the involved bonds breaks. These simulations 



are similar to those used in studies of lattice materials [ll|, ll3|, [l^ . We create a constant strain in the CRN using a 



constant displacement (A) and seed the system with an initial finite crack. Unlike the binary-alloy model for glass, 
where the crack always arrests under constant driving [l9l-[2l|, in our CRN model, the crack, except for the smallest 
drivings, always propagates through to the end of the sample, thereby yielding a steady-state crack for small and 
moderate driving (Fig. [IJa)). Only close to the threshold Griffith's displacement, Ac, does the crack occasionally 
fail to propagate. 

When we increase the driving displacement, micro-branching starts to appear (Fig. [Hb)). The simulation shown in 
Fig. [1] demonstrates an important feature of fracture in amorphous materials. Unlike lattice materials, in amorphous 
mater ialt he crack meanders off the midline. When the crack reaches a "tough zone" , where there is no bond in the 
preferred direction for breaking, i.e., back to the midline, the crack starts to bifurcate (Fig. [TJa) in the circle), but 
significant micro-branching will appear only if there is enough stored energy in the mesh, i.e. Fig. [IJb). There does 
not appear to be a sharp threshold for micro-branching; rather, the extent of micro-branching increases continuously 
with the driving. 

It turns out that the precise chosen crack path changes significantly with minor changes of our simulations, while 
the macroscopic features (the v{A) curve, the length of the micro-branches, etc.) remains very much the same, albeit 
with some noise. We took advantage of this, repeating the molecular dynamic simulations with the same driving 
conditions several times (typically 20), changing the Euler time-step each time by 3%, and averaging. Following the 
main cracky we measured the normalized velocity of the crack '^/cr as a function of the normalized driving displacement 
^/Ag. In Fig [21 we plot the '^(A) curve for two (amorphous) CRN's with different widths and for a lattice material. 
We can see that despite the fact that the two CRN's are completely different on the microscopic level, their v{A) 
curves are quite close to each other. They are also of the same magnitude as the v{A) curve of the corresponding 
honeycomb lattice material with the same force laws. In the lattice material, though, the critical velocity (the velocity 
when a steady-state crack stops propagating on the midline of the sample and starts to create micro-branching) is 
approximately 0.67ci^. In the CRN's there is no clearly defined threshold, and simulations exhibit micro-branching 
for A/ Ac well below the lattice critical velocity. This result is directly analogous to the experimental finding that 
the critical velocity of amorphous materials is much lower than that of ordered materials [6, 30]. 

As seen in Fig. 3, the CRN's behavior beyond the point of instability compares well with that of a propagating 
crack in experiments on a mode-I fracture in amorphous materials jH, [H, [3l|, as seen in Fig. [3l Fig. E[a) shows a 
photo of an experiment in PMMA. Micro-branching is seen clearly in addition to the main crack. Fig. [IJd) shows a 
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FIG. 2: (color online)The velocity of the crack, normalized to the Rayleigh velocity c^, as a function of driving displacement, 
A, normalized to the critical Griffith displacement for crack propagation, Ag, in two amorphous CRN's with different widths 
(82 and 162 atoms in the ^-direction, correspondingly) and in a honeycomb lattice with the same force laws. The dashed line 
represents the critical velocity of the lattice system. The amorphous systems have no well-defined threshold for micro-branching, 
and exhibit micro-branching for A/ Ac well below the critical velocity of the lattice system. In all cases, r] — 1.5. 



(b) 




FIG. 3: (color online) (a) An image of micro-branching in PMMA. The figure is taken from [3]. (b)-(c) Images of macro- 
branching in experiments. Image (b) is in Homalite-100 (glassy polymer) and taken from ^3^\ and (c) is in PMMA and taken 
from i. (d) An image of micro-branching fracture in our CRN simulation. 77 = 1.5, ^/ag = 3.1. (e) In higher driving (77 = 1.5, 
^/ag = 4), macro-branching begins to appear in our CRN simulation; both branches of the crack reach the end of the sample. 




parallel image from our simulation (of course, at a different length scale; our simulation is atomistic). Qualitatively, 
the images look quite similar. In the experiments, increasing the driving causes a bifurcation of the main crack to 
several (at least two) macro-branches traveling to the end of the sample, both in Homalite-100 (Fig. [3l^b)), and in 
PMMA (Fig. [3l^c)). Exactly the same feature is seen in our simulations (Fig. [3fe)). 

One of the interesting features revealed in the experiments is the purported increase in the amplitude of oscillations 
of the crack velocity with increasing driving (or increasing average velocity). In the steady-state regime, the v{t) curve 
is quite smooth. When the driving increases and micro-branching starts to appear, oscillations begin to appear in the 
electrical resistance of a conductive layer; attached to the sample [3', [HI, [6|. This was then interpreted as an oscillation 
of the crack velocity, based on a calibration of the resistance measurement, either by simulation of the conductive 
layer [si'-'si^, or experimentally by quasi-static crack measurements of a single straight crack [sj. Our simulations do 
not show an increase in oscillations of the velocity of the main crack as the crack velocity increases. Instead, we found 
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ir|=0.2, A/A^=2.7 
r|=0.15, A/A^=3.] 
:r|=0.15, A/A^=4 
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FIG. 4: (color online) (a)The rate of change of fracture surface, L, normalized to the Rayleigh wave speed, cr, as a function of 
the time for different drivings. For small driving, the curve is quite flat. At larger driving the curve starts to oscillate, due to 
the appearance of micro or macro-branching, (b) Correlation test to the crack's velocity oscillations. The oscillations of the 
electrical resistance behave like L{t), and not like the main crack velocity. Here r] — 1.5, A/Ag — 3.5. For both panels, the 
width is 162 atoms. 



CD 

o 

-1—' 

"co 

CD 

DC 
CO 

DC 



0.4 



0.2 



0.1 



• 










■ 


Narrow 
Wide 





























































2.6 



3.2 



3.4 



A/ 



3.6 



FIG. 5: (color online) The amplitude of the RMS fluctuations in the oscillations of the electrical resistance during the course of 
the simulation as a function of the normailized driving A/Ag- The "wide" system is 162 atoms wide, and the narrow one 82 
atoms wide. In both cases, 77 = 1.5. 



increasing oscillations in the rate of change of the overall length of the fracture surface, L{t) where L{t) = e Tii(t) 
and rii is one for a cracked bond. In Fig. Hl^a) we present three time histories of this new velociity. The black curve is 
for an essentially steady-state crack, the red curve corresponds to a crack with a small degree of micro-branching (Fig. 
[Sfd)) and the green curve is for a case when macro-branches were created (Fig. [Sl^e)). We can clearly see that the 
steady-state curve is much smoother than the other curves. In addition, the shape of the oscillations in the velocity 
curves are very similar to what was seen in the experiments [3l,l5l,[6|. 

It is straightforward to model the conductive layer by a grid of identical resistors, with the resistors removed 
whenever they cross a broken bond. Then for a specific crack configuration in our CRN we solve a Laplace-like 
equation for the voltage in the system, from which we can calculate the current density, and determine the sample's 
electrical resistance. A calculation of the resistance in our cracked CRN shows that it behaves like L{t)^ and not like 
the main crack's velocity (Fig. HJ^b)). Plotting the RMS amplitude of the oscillations of the electrical resistance (Fig. 
[5]) shows significant growth of the oscillations with the driving displacement, as in the experiments [3]. 

In conclusion, the main features of the micro-branching instability at large driving observed in the experiments 
were reproduced in our CRN-grid simulations both qualitatively (Fig. [3]) and semi-quantitatively (Fig. [5]). We 
should note that a fundamental question that our simulations do not resolve is the role of the third dimension in 
the experiments. The experiments on gels [35], show systematic changes with the thickness of the gel, indicating the 
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possible importance of such effects. We hope to address this issue after extending our model to three dimensions, 
which should in principle be straightforward, although computationally demanding. Also in the future, we plan to 
investigate the scaling properties of our CRN cracks as the system size is scaled up. 
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